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A Meta-Algorithm for Creating Fast Algorithms for Counting ON Cells 

in Odd-Rule Cellular Automata 


By Shalosh B. EKHAD, N. J. A. SLOANE, and Doron ZEILBERGER 


Abstract: By using the methods of Rowland and Zeilberger (2014), we develop a meta-algorithm 
that, given a polynomial (in one or more variables), and a prime p, produces a fast (logarithmic 
time) algorithm that takes a positive integer n and outputs the number of times each residue class 
modulo p appears as a coefficient when the polynomial is raised to the power n and the coefficients 
are read modulo p. When p = 2, this has applications to counting the ON cells in certain “Odd- 
Rule” cellular automata. (This article is accompanied by a Maple package, CAcount, as well as 
numerous examples of input and output files, all of which can be obtained from the web page for 
this article: 

http://WWW.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/CAcount.html). 

Preface 


The number of ON cells in the nth generation of an “Odd-Rule” cellular automaton is found by 
raising the defining polynomial (in which the number of variables is equal to the dimension of the 
ambient space) to the nth power, reading the coefficients modulo 2, and counting the remaining 
monomials—or equivalently, setting all the variables equal to 1 (see [SI] for a detailed discussion). 

The purpose of this article is to describe a meta-algorithm, inspired by a recent paper of Eric 
Rowland and Doron Zeilberger [RZ], that takes such a polynomial as input, and outputs a recurrence 
scheme that enables the fast (logarithmic time) computation of terms of the sequence giving the 
number of ON cells at time n. This provides an alternative, computer proof of Theorems 4 and 5 
of [SI]. 

A toy example 

Following the Gelfand Principle, let’s illustrate the method with a simple example that can be done 
by hand. We will later describe how this method can be ‘taught’ to a computer, which will then 
be able to do far more complicated cases, impossible for humans. 


Consider the sequence 

ai(n) := (1-|-X-|-x^)"" mod2 , 

X = 1 

(sequence A071053 in [OEIS]), and suppose we wish to compute ai(10^‘^‘^), or ai(n) for any very 
large n. 


Of course, direct computation is hopeless, even if we reduce modulo 2 at each step and use the 
repeated squaring trick that makes RSA possible (P” = (P"-/^)^ if n is even, P"" = pp^-^ if n is 
odd), since the polynomials, before we set x = 1, are far too big for our modest universe. What we 
will do is adapt this trick so that we can also make the substitution x = 1 at intermediate steps. 
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First let’s try to relate ai(2n) to ai(n), using the Freshman’s Dream identity P{xY = -P(x^) modp: 


ai(2n) = (1 + x + mod2 = ((1 + x + x^)^)” mod2 

X = 1 


= (1 + 2 = (1 + X + mod 2 

a:=l 


(replacing x^ by x). Hence 


ai(2n) = ai(n) 


{EvenCasel) 


(Recurrenceleven) 


Now we do the same thing for ai( 2 n + 1 ); 

ai( 2 n + l) = (1 + X + mod 2 = (1 + x + x^) ((1 + x + x^)^)"" mod 2 

X=1 X=1 

= (1 + X + x^) (1 + x^ + x^)"" mod 2 

X=1 

= (1 + x^) (1 + x^ + x^)’* mod2 + X (1 + x^ + x"^)"" mod2 . (OddCasel) 

X=1 X=1 

In the first term, once again, we can replace x^ by x, getting an uninvited guest, 02(n), say: 

a2(n) ;= (1 + x) (1 + X + x^)"" mod 2 

X=1 

As for the second term of Eq. (OddCasel), multiplying by x does not change anything, so this is 
equal to (1 + x^ + x"^)"" mod 2 , which, again replacing x^ by x, is our old friend ai(n). Hence 

X=1 

ai(2n + 1 ) = 02(71) + ai(n) . (Recur rencelodd) 

But this pair of recurrences is useless unless we can handle 02(0). So let’s try the same technique 
on it. A priori, this may force us to introduce terms 03(0), 04(71), etc., and lead us into an infinite 
regression, also known as a Ponzi scheme, but let’s hope for the best. 

Again we start with 02(277). Using the Freshman’s Dream, and the fact that multiplying a polyno¬ 
mial by X (or any other monomial) does not affect the result if we are going to read it modulo 2 
and set X = 1 , we have 

02(277) = ( 1 -|-x) ( 1 -|-X-|-x^)^"" mod 2 = ( 1 -|-x) • (( 1 -|-X-|-x^)^)” mod 2 

X=1 X=1 

= (1 -|- x) • (1 + x^ -|- x"^)” mod 2 = 1 • (1 + x^ -|- x"^)” mod 2 -|- x (1 + x^ -|- x"^)"" mod 2 

X = 1 X = 1 X = 1 

= 2 ( 1 -I-x^-I-x^)’* mod 2 = 2 ( 1 -|-x-|-x^)’* mod 2 = 204(77) . 

X=1 X=1 

Hence 

02 ( 277 ) = 201 ( 77 ) . (Recurrence2even) 


Now for 02(277 + 1). We have 

02(277 + 1) = (1 -|- x) • (1 -|- X -|- x^)^""^^ mod 2 

X=1 
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= ((1 + x) • (1 + X + x^)) • ((1 + X + mod2 

= (1 + 2x + 2x^ + x^) • (1 + x^ + ™od 2 
= (1 + x^) • (1 + x^ + x‘^)^ mod 2 


X = 1 


X = 1 


I X=-l 

J2 I _4\n . 


= 1 • (1 + + x^)'^ mod 2 + x"^ • (1 + x^ + x'^)^ mod 2 

x=l 


x=l 


= (1 + x"^ + x^)’^ mod 2 + (1 + x"^ + x^)’^ mod 2 

a:=l 


x=l 


= (1 + X + x^)” mod 2 + (1 + X + x^)"" mod 2 = 2 ai(n) . 

X=1 X=1 


Hence 


a 2 ( 2 n + 1) = 2ai(n) . {Recur rence2odd) 

So the uninvited guest, 02 (n), did not invite further guests, and now we have a super-fast way to 
compute ai(n) for large n, using the system 

ai( 2 n)=ai(n) , ai( 2 n-|- 1 ) = ai(n)-|-a 2 (n) ; 

a 2 ( 2 n) = 2ai(n) , a 2 ( 2 n-|-1) = 2ai(n) . (System) 


For certain “odd-rule” cellular automata, the sequence ai{n),n > 0 is completely determined by 
the subsequence hi{k) := ai( 2 ^ — l),fe > 0 [SI], and the bi{k), unlike the ai(n), often have simple 
generating functions, which we can derive (rigorously) by these methods. With ai(n) as defined 
above, let 

00 

fi{t) ■=^bi{k)t’^ 

fc =0 

be the generating function for bi{k), and similarly define 62 (fc) := ^ 2 ( 2 ^ — 1 ) and 

00 

f2{t) ■='^b2{k)t^ . 

fc =0 


From Eq. (System), we have 

b^(k) = b^(k-l) + b2(k-l) , b2(k) = 2b^(k - 1) , 


and since by direct computation, 61 ( 0 ) = 1 , 62 ( 0 ) = 2 , we arrive at a system of two linear equations 
for the unknowns /i (t) and /2 (t): 

{fl(t) = 'i-+tfl(t)+tf2(t) , f2(t) =2 + 2tfi(t)} , 


whose solution is 


hit) 


1 -F 2t 

(l + t)(l-2t) 


f2(t) 


2 

(l + t)(l-2t) 
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(A001045, A014113 in [OEIS]). But we really don’t care about / 2 (t), we just needed it in order to 
find so now we can safely discard it, and get the 

Theorem: 

f (f) = ^ 

(l + t)(l-2t) ■ 

The general case 

Fix once and for all a prime p and a polynomial P = P{xi,..., Xk) € Z[xi,..., Xk]- If A(xi,..., Xk) 
is any element of Z[a;i,..., Xk], we define the functional 

A(xi,..., Xk) —> A(a;i,..., Xk) modp (Reduce) 

Xi = l,...,Xk=l 

to mean “expand A(xi,..., Xk) as a sum of monomials, reduce the coefficients modulo p to one of 
the numbers {0,1,... ,p — 1} G Z, and finally set all the variables Xi equal to 1”. 

For any polynomial Q = Q(xi ,..., Xk) G Z[xi,..., Xk] whose degree in each of the variables is less 
than p, define 

aQ(n) := QP^ modp 

Xi = l,...,Xk = l 

For 0 < i < p, we have 

aQ{pn + i) = Q{xi,... ,Xk)Pixi,... ,Xk)^"'~^'' modp 

Xi=l,...,Xk = l 

= [Q{xi,.. .,Xk)P(xi ,... ,Xfc)*]P(xi,... ,Xfc)”'^modp 
= [Q(xi,... ,Xfc)P(xi,... ,Xfc)*](P(xi,... ,Xfc)P)”'modp 
= [Q(xi,...,Xfc)P(xi,...,Xfc)*]P(x^,...,x^)”'modp 

Now write 

Q(xi,...,Xfc)P(xi,...,Xfc)*modp = ^ x“" • • • (xf,..., x^) . 

(Here again “mod p” applies just to the coefficients, not the variables.) Hence 
aQ(np+i)= x“^ •••x“''i?(„^,...,„j^)(x?,...,x^)P(x?,...,x^)’"modp 

(Q:i,...,Q:fc)e{0,...,p-l}'= 

(Qi,...,afc)e{0,...,p-1}'= 
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Y1 • 

(ai,...,Q:fc)e{0,...,p-l}'= 

In other words for any Q{xi,... ,Xk) and each of the residue classes i,0<i<p—l,we can find a 
multiset of polynomials, let’s call it Si{Q), such that 

aqinp + i) = ^ aR{n) . 

ReSiiQ) 

We really only care about the case <5 = 1, but the algebra forces us to consider other Q’s, and they 
in turn force us to treat still other Q’s, and so on. However, by the pigeon-hole principle, this 
process must terminate, and we obtain a finite recurrence scheme, containing say m equations. 
Placing all the Q’s that appear into some arbitrary order, with Qi = 1, we get a (logarithmic-time) 

recurrence scheme: 

aj{np + i) = ^ ai{n) , 

for 1 < j < m, that enables the fast calculation of ai(n) for any n. 

Furthermore, by focusing only on i = p— 1, and defining Cj{k) := aj{p^ — 1), we have, for 1 < j < m, 

ci{k-l) . 

Define the generating functions 

OO 

fji*) ■= ^Cj{k)t’" (1 < j < m) . 

fc=0 

Standard manipulations of generating functions convert the above recurrences into a system of m 
linear equations for the m unknowns , fm(t)'■ 

fj{t) = Cj{0)+t ^ fi{t) , l<j<m , 

that can be solved, at least in principle, yielding rigorous explicit expressions for all the fj{t), 
and in particular for the one in which we are most interested. Note that this proves that 

the generating function, is always a rational function. If m is too large, and the system of 

equations cannot be solved, then one may try to use the recurrences to generate sufficiently many 
terms of the sequence ci{k), and then guess the rational function using for example the Maple 
packgage gfun [SaZ]. It may then be possible to justify that guess, a posteriori, by finding upper 
bounds on the degree of the generating function. 

Keeping track of the individual coefficients 

If instead of the functional Eq. (Reduce), one uses, for some formal variables si,..., Sp_i, 

> 
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one can modify the above arguments and keep track of the number of occurrences of each i (i = 

1,... ,p — 1) as coefficients in the expansion of P{xi, ..., Xk)^ modp. 

The Maple package CAcount 

Everything discussed above is implemented in the Maple package CAcount, which can be down¬ 
loaded from the web page for this article: 

http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/CAcount.html , where there 

are also many samples of input and output files that readers can use as templates for further com¬ 
putations. 

To see the list of the main procedures, type 
ezraO; , 

or to see the list of procedures that handle the more refined case, where one keeps track of the 
individual coefficients (only useful for p > 2), type 

ezraGO ; 

To get instructions on using a particular procedure, type 
ezra(ProcedureName); 

For example, procedure CAaut finds the recurrence ‘automaton’, and to get help with it, type 
ezra(CAaut); 

For our toy example, type 
CAaut([l+x+x**2,1], [x],2,2); 
which produces as output the pair 
[[[[1], [2, 1]], [[1, 1], [1, 1]]], [1, 2]] , 
where the first component, 

[[[1], [2, 1]], [[1, 1], [1, 1]]] , 

is Maple’s way of encoding the recurrence 

ai(2n) = ai(n) , ai(2n-|-l) = a 2 (n)-|-ai(n) ; a 2 ( 2 n) = ai(n)-|-ai(n) , a 2 ( 2 n-|-l) = ai(n)-|-ai(n) 


The second component 

[ 1 , 2 ] 
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is Maple’s way of encoding the initial conditions 


ai(l) = l , 02(1) = 2 . 

Procedure SeqF uses the scheme, once found, to compute as many terms as desired, while procedure 
ARLT (for anti-run-length-transform, see [SI]) computes the sparse subsequence in the places p* — 1. 
Procedure GFsP finds the proved generating function for that subsequence, and if the size of the 
system is too big, GFsG guesses it faster, and as we mentioned above, the guess can be justified a 
posteriori. 
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